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Abstract 

We have used Monte Carlo simulations, combined with finite-size scaling and two different real- 
space renormalization group approaches, to study a fully frustrated three-dimensional XY model on 
a simple cubic lattice. This model corresponds to a lattice of Josephson-coupled superconducting 
grains in an applied magnetic field H = (<3?o/a 2 )(l/2, 1/2, 1/2). We find that the model has a 
continuous phase transition with critical temperature T c = 0.681J//cb, where J is the XY coupling 
constant, and critical exponents a/v = 0.87 ± 0.01, v/v = 0.82 ± 0.01, and v = 0.72 ± 0.07, where 
a, v, and v describe the critical behavior of the specific heat, helicity modulus, and correlation 
length. We briefly compare our results with other studies of this model, and with a mean-field 
approximation . 

PACS numbers: 64.60.Fr, 74.81. Fa, 05.10.Ln, 64.60.Ak 
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I. INTRODUCTION 



The classical XY model has been widely used for decades as a model for phase transitions 
in materials with interacting spins. The dimensionality n of the spins (n = 2 for the XY 
model) is independent of the lattice dimensionality d, which can be either 2 or 3 for physically 
relevant systems. The d = 2 XY model undergoes the well-known Kosterlitz-Thouless (KT) 
transition [lj], characterized by an unbinding of vortex-antivortex pairs at the KT transition 
temperature Tkt- The KT transition is a continuous phase transition, but with unique 
critical properties 111]. The d = 3 XY model exhibits a more conventional phase transition 

u n 

with well-known critical exponents |2J, |3J. It is thought to describe many ferromagnetic 
materials with two-component spins. In addition, it describes phase transitions in which the 
"spins" actually represent the phases of a complex order parameter, such as the A transition 
in He 4 and superconductor-to-normal phase transition at zero applied magnetic field. 

The frustrated classical XY model, in either d — 2 or d — 3, has a much wider range 
of phase diagrams than does the unfrustrated case just mentioned. In this model, the 
coupling between spins is such that the ground state of the system cannot minimize all 
the bond energies simultaneously. Interest in this model was greatly increased when it 
was realized that this model described real systems, such as Josephson junction arrays in 
an applied magnetic field. The first demonstration that the fully frustrated XY model 
undergoes a continuous phase transition in d = 2 was given by Teitel and Jayaprakash j^J], 
using Monte Carlo (MC) techniques. This work was later extended to other values of the 
so-called frustration parameter / jjj, leading to an extensive literature on^two-dimensional 



(2D) frustrated XY model on various lattices and at different values of / 

The d = 3 frustrated XY model has also been studied extensively, in part because it is 
believed to describe flux line lattice (FLL) melting [s| under an applied magnetic field in 
high-T c superconductors. Hetzel et al. nj] used a uniformly frustrated 3D XY model on a 
stacked triangular lattice to study the melting of an unpinned Abrikosov lattice in a type-II 
superconductor. They showed convincingly that this melting transition is first order, rather 
than continuous — a prediction subsequently confirmed by experiment. Earlier work on a 



frustrated XY model on a simple cubic lattice, with magnetic fie 
lattice axes 



vy moaei on a simple cudic lattice, witn magnetic neia para 
Q found a continuous phase transition. Li and Teite, Q Q 



d parallel to one of the 
used a uniformly 



frustrated XY model similar to that of Ref. [10] to calculate the properties of the vortex 
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line liquid which appears above the melting temperature, and to investigate the possibility 
of a further phase transition between an entangled and a disentangled vortex line liquid 
above the FLL melting transition. Chen and Teitel 1J| later extended this work to the more 
realistic case of uniaxially anisotropic couplings; they suggested that there was another phase 
transition temperature T cz above the FLL melting temperature T m , where superconducting 
coherence parallel to the applied magnetic field would vanish. Chin et al. jl^ have, however, 
suggested that this apparent existence of a transition to a disentangled vortex liquid is a 
result of finite system sizes and simulation times. More recently, the uniformly frustrated 3D 
XY model was studied by Monte Carlo methods on a simple cubic lattice with a magnetic 

n 

field parallel to the [111] direction [161 ] . For this choice of field direction, the simple cubic 
lattice behaves as a stack of 2D triangular lattices with ABC ABC ■ ■ ■ stacking. It was 
found that this system, similar to that studied in Ref. Q], exhibits a clear first-order FLL 



melting transition. Nguyen and Sudb0 [17| considered a uniformly frustrated anisotropic 
Villain model (an approximation to the uniformly frustrated XY model) to study the phase 
diagram of a uniaxially anisotropic high-T c superconductor as a function of the applied 
magnetic field and temperature. They found two phase transitions: the lower-temperature 
one is FLL melting, while that at higher temperature involves the destruction of the phase 
coherence in the direction of the applied magnetic field. 

In this paper, we study phase transitions in the fully frustrated XY model in d — 3, using 
primarily MC simulations. This model is of interest, in part, because it may be relevant to 
FLL melting in the high-T c materials. In addition, because of great advances in microfab- 
rication techniques, it is now possible to make microscale or nanoscale arrays of Josephson 
superconducting grains in three dimensions. Such an array should, in a suitable applied 
magnetic field, be describable by a frustrated XY model, at least to a first approximation. 
It has also recently been suggested that a fully frustrated XY model might also be realized 
by an assembly of cold atoms on a suitably constructed optical lattice jlja]. 

The fully frustrated XY model is characterized by the frustration vector f = 
(1/2,1/2,1/2), as further defined below. Such a model has been previously studied by 
Diep et al. [19], who found that, in contrast to the f = (1/3, 1/3, 1/3) case studied in Ref. 
1^ | , there was a continuous phase transition. We extend the work of Ref. [3] by calculating 
the critical behavior of the helicity modulus, equivalent to the spin-wave stiffness constant 
(or to the superfluid density in a superconductor). We also carry out a more extensive finite- 
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size scaling analysis than done by those workers, thus obtaining more accurate information 
about the critical behavior. Our results do, however, confirm that the phase transition is 
continuous, not first order. 

The remainder of this paper is organized as follows. In Sec. |H1 we present the formalism 
for calculating the thermodynamic properties and the critical exponents of our model. In 
Sec. IIII1 we give the results of our Monte Carlo calculations, finite-size scaling methods, and 
renormalization group methods. Section ITVl presents a summary and discussions. 

II. FORMALISM 

We now describe the model Hamiltonian on which our calculations are carried out. For 
convenience, we present this Hamiltonian as it applies to a simple cubic lattice of supercon- 
ducting grains in the presence of an applied external magnetic field H, though the model 
is not limited to this application, of course. Then the frustrated XY model in d = 3 is 
described by the Hamiltonian 

7-C = Jij cos(0i - - A {j ). (1) 

(ij) 

Here (pi is the phase of the superconducting order parameter on the ith site, Aij = 
(27r/$o) // A ' dl, $o = hc/2e is the flux quantum, A is the vector potential, is the 
coupling constant between the ith site and jth site, and (ij) represents a sum over all dis- 
tinct pairs of nearest-neighbor sites on a simple cubic lattice. We assume a constant coupling 
between each nearest-neighbor pair of sites, so = J and J > 0; we neglect the possible 
dependence of J on the ap plie d magnetic field H and the temperature T. We also assume 
weak screening as in Ref. [16]. Thus, the local magnetic field B is approximated by the 
applied magnetic field H. 

The x component of the frustration vector f = (f x , f y , f z ) is defined by 

J2 {x) A 3 =2nf x , (2) 

p 

where the sum is taken along the sides of a plaquette on the yz plane of the lattice; analogous 
definitions hold for f y and f z . If the simple cubic lattice has lattice constant a, fa = -Bja 2 /$o 
represents the flux through a single plaquette perpendicular to the ith axis, in units of one 
flux quantum. 
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One can also define the vortex number (or vorticity) of each square plaquette. For exam- 
ple, the vorticity of a plaquette lying in the yz plane is defined by 

n x = f x + {X) - & ~ A ^ ' ( 3 ) 

p 

where the sum runs counterclockwise around the perimeter of the plaquette, viewed from 
the positive x direction, and the phase differences are chosen so that < 4>i < 2n. n y and 
n z are defined analogously. 

In order to study possible phase transitions within this model, we have calculated several 
quantities. One of these is the specific heat Cy per site, given by 

where N is the total number of sites in the lattice, Ti is the Hamiltonian in Eq. (JI}, and 
(• • • ) denotes an average within the canonical ensemble. A first-order phase transition is 
generally indicated by a 5-function-like anomaly in Cy, while a continuous phase transition 
is signaled by lattice-size-dependent divergence in the Cy. 

To study the vortex lattice melting, we calculate a suitable vorticity density- density 
correlation function. Specifically, we first introduce the Fourier transform of the vorticity 
density ra*(k) = J^ R nj(R) exp(zk ■ R), where n»(R) represents the vorticity of a plaquette 
centered at R and oriented perpendicular to the zth axis. We then calculate some of the 
correlation functions 

R,R' 

where the sum is carried out over all R and R' such that R' — R = r. In our actual 
simulations, we use periodic boundary conditions in all three directions. In practice, it 
is more convenient to compute the Fourier transform of <?y(r). This Fourier transform is 
known as the vortex structure factor, and is given by SV,'(k) = g i:; -(r)e ik ' r . Because of the 
periodic boundary conditions, SV,(k) is defined only for k = [2ir / (N x a)](rni, ra 2 , rn 3 ), where 
rrii are integers, and the computational cell is assumed to contain iV = iVj sites. We obtain 
gij(v) by first computing 5V; (k) and then Fourier-transforming back into real space. 

The phase transition in this model is best characterized by the helicity modulus tensor 
7a/3 jscj- la/3 measures the stiffness of the phase against an external twist. It is defined as 
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the second derivative of the free energy with respect to an infinitesimal phase twist |2l| 



la(3 



d 2 F 



(6) 



5=0 



N d5 a d5 p 

7 Q/ 3 is conveniently calculated by adding a fictitious vector potential A' to the Hamiltonian 
in Eq. (JTJ), as 



d 2 F 



N dA'dA'- 

1 J 



(7) 



A'=0 



in the presence of periodic boundary conditions on the phases. The resulting diagonal 
components are readily evaluated, with the result 



7o 



+ 




Nk B T 



(8) 



Here e^ is the unit vector from the ith site to the jth site, and e a is the unit vector in the a 
direction. Since the helicity modulus 7 aa is proportional to the a component of the superfluid 
density, it follows that when ^ aa > 0, there is nonzero phase coherence in the a direction, 
and when ^ aa — > 0, the phase coherence is lost. Therefore, the superconductor-to-normal 
phase transition occurs at the temperature at which ^ aa — > 0. 

The correlation time r is a measure of how long it takes for the system to lose its memory 
of its previous state, r can be obtained by calculating a time-displaced autocorrelation 
function |22|. The time-displaced autocorrelation function xe(^) of the energy, for example, 
is defined by 



XE (t) = j[E(t')-{E))[E(t' + t)-{E))dt> 
= j[E(t')E(t' + t)-(E) 2 }dt', 



(9) 



where t and t' are two different Monte Carlo "times." Since Monte Carlo simulations involve 
fictitious dynamics, these times are not related in any obvious way to a physical time. 
However, the length of the Monte Carlo time does give a measure of the persistence of 
various MC states in time. 
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Equation Q can be expressed more conveniently when we have a set of measurements 
of the energy E(t) from t = (after equilibration) up to some maximum time t max . In this 
case, Eq. Q becomes 



^max t 



XE(t) = —— £ E(t')E(t' + t) 



^max t ^max t 



I "max " "max " 

-— - £ s(0 £ £(;' + ;). (io) 



(tmax t) 2 t , =Q t/=Q 

Since the autocorrelation function is expected to fall off exponentially at long times as 

XiK*)~e-*/ r , (11) 
the correlation time r can be calculated from 

x E (t) 



dt = / e-* /r dt = r. (12) 



/o Xs(0) 7o 

This expression for r is also called integrated correlation time. MC "measurements" will be 
statistically independent only if they are separated by intervals of ~ 2r or more. In addition, 
the correlation time r is related to the dynamic exponent z by 

r ~ i z ~ L z , (13) 

where £ is the correlation length (equal to L at the T C (L)). The dynamic exponent measures 
the extent of the critical slowing down. The smaller z is, the more accurate are the numerical 
measurements. 

If there is a continuous phase transition, various quantities should exhibit singular be- 
havior near the transition temperature T c . If we define a reduced temperature by 

T — T 

J- c 

then in the thermodynamic limit the correlation length £, the specific heat per site Cy, and 
the helicity modulus 7 near T c are expected to vary as 

e ~ i*r, (15) 

Cy ~ (16) 

7 ~ \t\ v , (17) 
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where v, a, and v are critical exponents. If T c is known, the simulation data can be fitted 
to the expected asymptotic form to obtain values of the critical exponents. But since T c 
is typically not known in advance, it is usually more accurate to use a different method. 
One such method is finite-size scaling, which extracts values for the critical exponents by 
investigating how measurements depend on the size L of the system [23]. This procedure is 
carried out by expressing a quantity of interest in terms of the correlation length and then 
introducing a new dimensionless function, known as a scaling function. For example, Cy 
and 7 are expressed as 

C v (L,t) = L a ' u C v (L l / u t) } (18) 
7 (L,t) = L-^(L^t), (19) 

where L = N x a is the linear system dimension. Since the scaling functions Cy and 7 should 
depend on a single variable, we can make all the data for each system size L fall on the 
same curve by appropriately adjusting the values of the critical exponents and T c . When 
this happens, we assume that we have the correct values for these quantities. 

Another method is the phenomenological renormalization group (PRG) HQ |2J. To 



calculate the critical behavior of 7, for example, we consider two different system sizes L 
and V and introduce the ratio 

W> s 3£| (20) 

From Eq. (fT9"j) this ratio becomes P 7 = (L/L')- v / v when t = 0. Therefore, if one plots two 
different curves of P-y{L, L', t) versus t with the same ratio of L/L', the temperature at which 
they intersect is T c , and the value of P 7 at that temperature yields the ratio vjv. 

The critical temperature T c , or its inverse value K c = J / (/cbT c ), and the critical exponent 
v can also be determined from Binder's fourth-order cumulant Ul 0] defined by 



Ul = 1- vJ^, (21) 



where s = (1/N)y (^2f cos0j) 2 + (Xlf sin 0^) 2 . The scaling form for the Binder's cumulant 
is Ul = U[L x l v t\ without any prefactor. Hence, Ul can be Taylor expanded about T c as 

U L = Ua + UiL 1 /" 

= U + UiLV v (l-jf)+.--. (22) 



If we plot Ul for several values of L as a function of temperature or its inverse value, they 
will intersect at the critical temperature T c . To obtain the exponent u, we can calculate 
{dU L /dK) K=Kc , where K = J/{k B T). From Eq. ((22), we find that 

dU L 

k, K c 



dK 



(23) 



Hence, the ratio of the two slopes for different values of L gives v. 



III. MONTE CARLO CALCULATIONS 



To carry out our Monte Carlo calculations, we used the standard Metropolis algorithm 
with periodic boundary conditions in all three directions. We started with a random phase 
configuration at temperature T = l.OJ/ks, then cooled down to T = OAJ/ks in steps of 
0.01J//cb, except near T c , where we decreased the temperature in steps of 0.005 J/ k-g. At 
each T, we took 50000 Monte Carlo steps per site through the entire lattice for equilibration, 
and then calculated the expectation values of the quantities of interest by averaging over an 
additional 20000 MC steps. We also used the final configuration of the previous T as the 
starting one of the current T. Near T c , the system undergoes critical slowing down; so we 
increased the number of MC steps up to 5 x 10 5 for equilibration and 2 x 10 5 for averaging. 

Instead of considering continuous angles between and 27r for the phases 0, of the order 
parameter on each site, we used the 360-state clock model, which allows angles of 0°, 1°, 
2°,..., 359°. This simplification should have no effect on our results, since it is known that 
there is no distinction between the continuum and the discrete results for the n-state clock 
model when n > 20 [271 ]. 

To calculate the phase factors we use the gauge A = ($ /a 2 )(f y zx + f z xy + f x yz), 

(z) 

where the frustration is f = f x x + f y y + f z z. Thus, for example, the phase factors A\-' 
arising from the field component parallel to z all vanish except for bonds in the y direction; 
for these bonds, and x = na, J±f) is given by 

A$ = 2imf z . (24) 

The phase factors and A^f are given by analogous expressions. For f = (1/2, 1/2, 1/2) 
and periodic boundary conditions in all three directions, all the phase factors are equal to 
either or tt, and thus all the couplings are of the form ±Jcos(0j — 0,). This choice of 
gauge automatically satisfies the condition given in Eq. (0). 
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We turn now to our numerical results. The values of the correlation time r for several 
lattice sizes are given in Table I, at the estimated transition temperature T C (L) for the fully 
frustrated XY model [f = (1/2, 1/2, 1/2)] of size L = N x a. r clearly increases monotonically 
with increasing lattice size. From the fits to data in Fig. ^ we get the dynamic exponent 
z = 2.23 ± 0.14 in our system. 

The internal energy per site U = E/N, expressed in units of the coupling constant J, is 
shown in Fig. |21 It is quite size- dependent for L < 10 but quickly converges for L > 10. The 
sharp drop in U near T = 0.7J//cb suggests but obviously does not prove the occurrence of 
a phase transition near that temperature. 

The calculated specific heat per site Cy is shown in Fig. OJ There is a clear peak near 
T ~ 0.7J/A;b which becomes sharper with increasing L, suggesting a continuous phase 
transition. This behavior is similar to that of the 2D fully frustrated XY model as in Ref. 
} and that of the ordered simple cubic lattice with f = (0, 0, 1/2), as in Ref. |ll|. As n 
those two examples, the finite magnitude of Cy at its peak is a result of the finite system 
size; otherwise, Cy would diverge at T c in the infinite system. 

Next, we turn to the behavior of the helicity modulus tensor. Since ^ xx = ^ yy = ^ zz for 
this isotropic system, we calculated the average 7 = (j xx + jyy + 7 zz )/3 in order to improve 
the statistics. j(T) increases with decreasing T, and shows a fairly clear drop to near zero 
near T = 0.7J//cb when L > 10, as can be seen in Fig. HJ For L = 10 or smaller, j(T) shows 
a broad transition region. A more accurate value for T c will be given further below. 

To estimate the statistical errors in our calculation of j(T), we used the jackknife method 
[3] , which is carried out as follows. For each T, we made n independent numerical measure- 
ments of 7(T), with measurements separated by at least two correlation times. From these 
n measurements, we calculated a value j(T) for the helicity modulus. Next, we removed 
the first measurement from this set of n measurements to calculate the helicity modulus 71 
with the remaining n — 1 measurements. To calculate 72 we restored the first measurement 
to the set and removed the second measurement, and so on. Thus, 7» is calculated with the 
ith measurement removed from the set. The error estimate for 7 is given by 



The error bars from this method are shown in Fig. but they are smaller than the symbol 
sizes. 




(25) 
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In Fig. we show the evolution of the internal energy per site U as a function of MC 
time at the transition temperature T c = 0.681J//cb for a lattice size 12 x 12 x 12. The three 
intensity plots for each a and b show the density-density correlation functions g zz (x, y, L z /2), 
9xx(L x /2, y, z), and g yy (x, L y /2, z) (denoted by xy, yz, and zx, respectively) of the vortices 
at two different times as indicated in the energy evolution curve. Each intensity plot is an 
average of one correlation time r = 148 MC steps per site through the entire lattice. The 
two times correspond to energies slightly below and slightly above the average value and 
they are separated by a very large MC time. In contrast to the results of Ref. ^(|, there 
is no indication of a vortex lattice phase in these plots; instead, both seem to show vortex 
liquid phases (although there are partial latticelike formations, especially on the zx plane 
in window a). This behavior suggests (as does the diverging specific heat peak in Fig. |3J) 
that the phase transition at f = (1/2, 1/2, 1/2) is continuous, rather than first order as at 
f = (1/3,1/3,1/3). However, at very low temperatures (T = 0.01J/k B ), the correlation 
functions show a clearer indication of an ordered phase, as shown in Fig. Specifically, 
we see evidence of a checkerboard lattice in all three directions, but especially in the yz 
windows. 

This ordering is clearer if we look at the vortex structure factor S zz (k) of the lattice. 
Figure [7| shows this structure factor for k parallel to the magnetic field, at several tempera- 
tures, including the two shown in Figs. El and El Here kn = |k| and k = [2n / (N x a)](m, m, m), 
where m — 0, 1,..., N x — 1. At low temperatures, the vortex structure factors at k\\ = V3ir/a 
are nonzero, but they do not increase monotonically with decreasing temperature. Instead, 
they increase down to T = 0.30 J/ k B , then decrease down to T = 0.01J/k B . We believe that 
this behavior may represent some kind of "polycrystalline" domain structure of the vortex 
at low temperatures. As shown below, the system has several eigenmodes, and the exact 
admixture of these eigenmodes may change as the temperature varies. The insets to Fig. 
[7| are the vortex density-density correlation functions g zz (x,y, L z /2) at two different tem- 
peratures T = 0.30J//cb and T = 0.681 J/ A;b- We can see a clear checkerboard pattern at 
T = 0.30 J/ &B- The vortex structure factor appears to go to zero near the same temperature 
as T c where the helicity modulus vanishes, although we did not collect enough numerical 
data to determine the temperature dependence of the structure factor peak near T c . This 
point is discussed further below. 

In Fig. |SJ we show the probability distribution P(U) of the internal energy per site U for 
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several lattice sizes at T c . For each lattice size, we see only a single peak, which sharpens 
with increasing lattice size. This result is also in contrast to that of Ref. [lfj], where the single 
peak splits into two peaks as the lattice size increases. This persistent single-peak behavior 
suggests a continuous phase transition, in contrast to the first-order phase transition found 
in Ref. Q- 

The data shown in Fig. |3] can be used in a standard way to obtain an estimate of the 
critical exponent ratio at/ 'v. From Eqs. (JT5j) and ([16)1 . the maximum value of Cy for a lattice 
of edge L is 

In Fig. El we plot logC™ ax versus logL at T = T c . The data is well described by a straight 
line. The slope of this fitted straight line gives a/v — 0.87 ± 0.01. The present result is in 
contrast to the linear relation between Cy and logL, corresponding to a = (logarithmic 
divergence) seen in the fully frustrated 2D XY model as in Ref. ^J. 

As we did for the specific heat Cy, we see from Eqs. ()15j) and (|T7jl . or just from Eq. ()19|). 
that the helicity modulus 7 depends on the lattice size L as 
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~ L- v/u . (27) 



In Fig. E3 we plot log 7 versus logL at two different temperatures T = 0.681 J/k B and 
T = O.Q82J/k B . From the fits to the data, we get v/u = 0.79 ± 0.02 at T = 0.681 J/k B and 
vjv = 0.83 ± 0.02 at T = 0.682 J/k B . The error bars from the jackknife method are also 
shown, but they are smaller than the symbol sizes. The value of vjv is very sensitive to 
small changes in the assumed T c . So we need a more reliable method to calculate vjv. In 
Fig. HH we show the results of a PRG study of the helicity modulus j(L,t). As described 
near Eq. (J20J) , we plot P 7 (L,L',t) for the fixed ratio L/L' = 2, but for three different pairs 
of sizes L and V . From the intersection of these three -P 7 (L, L', t) curves as a function of t, 
we find T c = 0.681 ± 0.001 J/k B and v/u = 0.82 ± 0.01. 

The previous two methods give only the ratios of two critical exponents. In order to find 
all three critical exponents, we need a third method to extract v. For this, we use another 
renormalization group transformation We first calculate U(T) for al6xl6xl6 lattice. 
Next, we reduce each linear lattice dimension by a factor of 2, by combining a cube of eight 
adjacent sites into one; this is known as a blocking scheme of the real-space renormalization 
group method. This factor 2 corresponds to the scaling factor b. The phase of the merged 
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eight sites is assigned by a simple additive rule — that is, it is taken as 

b = arctan [ Sm - ) , (28) 
\Ei cos0,y 

where — 7r < 0b < tt- We add 27r when 0b < so that 0b satisfies < 0b < 27r. The 
correlation length £' of the blocked lattice is also reduced by the scaling factor 6, i.e. 

£' = | (29) 

when £' is measured in terms of the lattice constant. But since the correlation length should 
diverge at T c , it follows that 

£(T C )=£'(T C ). (30) 

To use this scheme to obtain u, we also calculate the internal energy per site U' of the 
rescaled 8x8x8 lattice as a function of temperature T. This rescaled lattice should behave 
very similarly to an 8 x 8 x 8 lattice at a different temperature T'\ so 

U\T) = U{T'). (31) 

If we rewrite T' in terms of T, we obtain 

T' = U- 1 [U'(T)}, (32) 

where U~ l is the functional inverse of the function U. This relation between T and T' is 
the desired renormalization group transformation. As in Eq. (|15jh the rescaled correlation 
length £' can be expressed in terms of the reduced temperature t' as 

~ \t'\- v . (33) 

From Eqs. (fTT)|l . ([25)1 . and (f33)t . we obtain 

^)"" = 6. (34) 

Therefore, knowledge of the renormalization group transformation taking the system from 
T to T' gives the value of v. 

Since Eqs. (}2"9"|) and (}3*3*j) are the asymptotic forms near T c , we can linearize the renor- 
malization group transformation using a Taylor series expansion near T c to obtain 

dT' 



T'-T C =(T-T C/ dj . 
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(35) 



or equivalently 

a dr 



t dT 

Inserting this result into Eq. (J34j) yields 

log b 



(36) 



log 



dT \T C 

In Fig. we show simulation results for the rescaled 8x8x8 lattice and for a sep- 
arate unrescaled lattice of the same size. If we fit the two sets of points near the critical 
temperature T c to two straight lines, the ratio of the slopes yields the value of (dT'/dT) Tc , 
according to Eq. (|36j) . We can then get v from Eq. (|37p . From the plot shown in Fig. we 
obtain v = 0.72 ± 0.07. Since we already know ajv and v/u, we can use this value of v to 
obtain a = 0.63 ± 0.07 and v = 0.59 ± 0.07. In practice, the plot for the rescaled data has 
considerable uncertainty in the slope; so v cannot be determined with great accuracy using 
this method, at least using simulations of the size we have carried out here. 

All of our numerical data appears to be consistent with a single phase transition at 

T c = 0.681J/A;b- By contrast, there are two separate phase transitions in two-dimensional 

fully frustrated XY model on a square or triangular lattice, as discussed in Refs. j30| and jjsi |. 

One of these is a KT transition while the other is in the Ising universality class. Although 

the two transition temperatures are close to each other, the Ising transition temperature is 

slightly higher than that of the KT transition. Reference 3Jj explains the sequence of these 

phase transitions in terms of the loss of phase coupling across a domain wall. Although 

our results are consistent with a single phase transition, we do not have sufficient data to 

rule out the possibility that there could be two separate phase transitions in our model, as 

further discussed below. 

To gain further insight into our numerical results, we have also used the mean-field 

I I 

approximation developed by Shih and Stroud |32| to find the relevant eigenmodes near the 
phase transition at T c . Near T c , the linearized mean-field equations are 

m ~ § E •V ''/, = °> ( 38 ) 

j 

where 77, = (e 1 ^) and /3 = l/k B T. If we assume = J and use T' = ksT/J, Eq. (|3"%|) 
becomes 

^-^E e ^' = - ( 39 ) 

3 
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In the mean-field approximation, T c is the highest value of T such that Eq. (}3"9"j) has a 
nontrivial solution. With our frustration f = (1/2,1/2,1/2), the Hamiltonian is periodic 
with a unit cell of 2 x 2 x 2, which implies that Eq. (J39|) is a set of eight coupled homoge- 
neous linear equations. T c is given by the condition that the determinant of the matrix of 
coefficients should vanish. This requirement leads to two values of X" = ±v^3/2 ~ 0.866; 
each is fourfold degenerate. Thus, the transition temperature in the mean-field approx- 
imation is T C MF = 0.866 J/k-B- The four degenerate eigenmodes corresponding to this 
eigenvalue are (0, -1, 1, \/3, 0, 0, 0, 1), (1, 0, ->/3, -1, 0, 0, 1, 0), (-1, V3, 0, -1, 0, 1, 0, 0), and 
(— y/3, 1, 1, 0, 1, 0, 0, 0); these modes are shown in Fig. ED Of course, any linear combination 
of these modes would also be an eigenmode. Since the four eigenmodes are degenerate, they 
must be related by some discrete symmetry operations of the lattice. Thus, we might expect 
that there are several degenerate ground states. Indeed, our simulations at very low tem- 
peratures (T = 0.01J//cb) do suggest that the system can readily fluctuate among several 
different states at such temperatures. 



IV. DISCUSSION 



We have investigated the phase transition in a fully frustrated 3D XY model on a simple 
cubic lattice, corresponding to an applied magnetic field H = (Q /a 2 )(l/2, 1/2, 1/2). In 
contrast to the case f = (1/3, 1/3, 1/3) as in Ref. we see a continuous phase transition. 
We also extract the critical exponent ratios ajv = 0.87± 0.01, vjv = 0.82 ± 0.01, and, with 
less accuracy, v itself, using a variety of numerical techniques. We get a = 0.63 ± 0.07, 
v = 0.59 ± 0.07, and v = 0.72 ± 0.07. 

It is of some interest to compare these values with those of other models. In the isotropic 
unfrustrated) 3D XY model, v ~ 0.66 - 0.67 Q, E^li 33 34 35, 36, 37, 38,0. « ~ 
3, H Q (a = -0.017 in Ref. Q), and v ^ J&U BF 

l^ol l reported that a 



. For the anisotropic but 
- -0.007. In the weakly 



unfrustrated XY model in d — 3, Ref. 
frustrated XY model in d = 3 (f = (0,0,/), with / < 1/12), it has been reported that 
v pa 1.5 (4^. Reference found that v = 2.2 ± 0.4 in a random- coupling 3D XY model 
with free boundary conditions and f = (0, 0, l/27r), while Ref. found that v — 1.1 ± 0.2 
in a random-coupling 3D XY model with periodic boundary conditions and f = (0,0, 1/4). 
Our value of v satisfies the usual trend that v becomes larger than 0.66 — 0.67 when the 
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magnetic field is nonzero. However, our results appear not to satisfy the two hyperscaling 
laws a = 2 — dv and v — (d — 2)v. Possibly the reason is that at this phase transition, two 
order parameters go to zero: the helicity modulus, and a discrete, Ising-like order parameter 
related to the amplitude of the structure factor at a characteristic wave vector. 

The T c we obtain from the PRG method for 7 agrees very well with that found from 
the RGT method, and quite well also with that estimated from Cy(T). It is, however, 
about 20% lower than that obtained by the mean-field approximation. This deviation is 
not a surprise, because the Monte Carlo T c has been shown to be lower than that found by 
mean-field theory in other frustrated XY systems 11|. 

In order to avoid the problem of critical slowing down, instead of just increasing the 
number of Monte Carlo steps, one might attempt to use the cluster flipping algorithm, 
as first proposed by Swendsen and Wang |0, and later by Wolff 0], and many others 



45 



46 



47 



48. 



49 



501 ]. This algorithm is very successful in reducing critical slowing down 



for unfrustrated models such as the 0(n) a model and the Potts model. But when the 
system is frustrated, especially fully frustrated, a mere extension of the cluster algorithm 
does not reduce the critical slowing down, and may even increase the correlation time jf|. 
The reason is that the cluster percolation temperature T p is much larger than T c in the 
frustrated system, rendering the cluster flipping trivial at T c because the percolating cluster 
takes up almost the entire system. Therefore, the cluster should be generated in such a way 
that T p is closer to T c [fj, 0]. Critical slowing down might also be reduced if the standard 
Metropolis algorithm is combined with the cluster algorithm, resulting in a hybrid algorithm 



vietropons 
3,0,0- 



The renormalization group transformation method is based on certain assumptions which 
may lead to systematic errors 29]. Specifically, it assumes that the blocked system has a 
typical phase configuration of another 3D XY model on a simple cubic lattice with V = L/b, 
i.e., that they appear with the correct Boltzmann probabilities as the original states. This 
assumption is not exactly correct, and contributes to some systematic error, which cannot 
be easily estimated. In the present paper, we have tried to estimate these errors using the 
jackknife method. The renormalization group method is also affected by finite-size effects, 
because the original system has a different size than does the rescaled system. To optimize 
the benefits of this method, therefore, we should ideally run simulations on as large a system 
as possible. We should also do an extra simulation on a system whose size is the same as 
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that of the rescaled system, and compare the two results. 

We also comment briefly on the difference between our work and that of Diep et al. 
These authors did not compute the helicity modulus, nor did they comment on the 
connection between the model and a superconducting array in a magnetic field. In addition, 
because of the several techniques described above, we are able to get more accurate values 
of the ratio a/u, as well as values for v/u and of u itself. 

Finally, we briefly discuss the possibility that the phase transition at T c might actually be 
two separate phase transitions. In the 2D fully frustrated XY model, as already mentioned, 
there are indeed two separate phase transitions: a KT transition at a lower temperature, 



followed by an Ising-like transition at a slightly higher temperature |30|, |3jj between a state 
in which the vortices are ordered in a checkerboard pattern and a disordered vortex state. 
In the present case, the transition at T c also involves the nearly simultaneous disappearance 
of 3D XY order (signaled by the vanishing of the helicity modulus) and a discrete order 
parameter (indicated by the vanishing of the vortex structure factor). Once again, in 3D, 
the discrete order is characterized by a checkerboard configuration of the vortex state below 
T c (see Figs. EHZI), although in this case the discrete order parameter is described by four 
degenerate modes, as determined by the mean-field solution, rather than two as in the 2D 
fully frustrated XY model jo^ . 

Although the discrete and XY order appear to vanish at the same temperature in 3D, and 
there is no evidence of two separate phase transitions, we believe that our numerical results 
are not sufficient to conclusively rule out two separate phase transitions. In particular, we 
have not carried out careful numerical studies of the critical behavior of the discrete order 
parameter. It would be of interest to carry out further numerical studies, especially of the 
discrete order parameter, to answer this question definitively. 

To summarize, we have studied the phase transition in the fully frustrated XY model, 
using Monte Carlo simulations in conjunction with two types of real-space renormalization 
group approaches. We find, in agreement with previous work, that the phase transition is 
continuous, and we obtain accurate values of the critical exponents a/u and v/u, and a 
slightly less accurate value for u itself. The phase transition could, in principle, be probed 
experimentally in a suitable three-dimensional lattice of coupled superconducting grains, 
and possibly also in an assembly of cold atoms in an optical lattice. 
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TABLE 



L 


t(*mc) 


4 


14 


6 


29 


8 


53 


10 


94 


12 


148 


14 


164 


16 


352 



TABLE I: Calculated correlation time r for several lattice sizes L, evaluated at T C (L). r is 
measured in the units of MC steps per site. 
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FIGURES 




FIG. 1: Log-log plot of the data in Table I, corresponding to a fit of r to the function r = AL Z . 
The slope of the fitting line to the data yields the dynamic exponent z = 2.23 ± 0.14. 
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FIG. 2: Internal energy per site as a function of temperature T for several lattice sizes, as indicated. 
Note that T decreases with increasing distance along the horizontal axis. 
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FIG. 3: The specific heat per site Cy as a function of temperature for several lattice sizes. The 
lines are cubic spline fits to the data. 
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FIG. 4: The averaged helicity modulus 7 = (j xx + ^ yy + 7zz)/3 as a function of temperature for 
several lattice sizes. 
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FIG. 5: Internal energy per site for a lattice of size 12 x 12 x 12 as a function of MC time, at the 
transition temperature T c = 0.681J/£;b- The intensity plots a and b represent the vortex density- 
density correlation functions g zz {x, y, L z /2), g xx (L x /2, y, z), and g yy (x, L y /2, z) denoted by xy, yz, 
and zx, respectively, at the times a and b in the energy evolution curve. A lighter color represents 
a larger value of the correlation function. 
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FIG. 7: Vortex structure factor S zz (k\\) parallel to the applied magnetic field at several tempera- 
tures for a lattice of size 12 x 12 x 12 and periodic boundary conditions. The insets are the vortex 
density-density correlation functions g zz (x,y, L z /2) at two different temperatures T = 0.30J/fcB 
(top) and T = 0.681J/£;b- A lighter color represents a larger value of the correlation function. 
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FIG. 8: MC probability distribution P(U) for the internal energy per site U at the transition 
temperature T c = 0.681J//cb for several lattice sizes. 
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FIG. 9: Log-log plot of the maximum height Cy of the specific heat per site versus linear size L 
of the lattice. The points are MC data; the full curve is the best-fit line. The slope of the fit line 
gives a/v = 0.87 ± 0.01. 
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FIG. 10: Log-log plot of the helicity modulus 7 versus linear size L of the lattice. The points 
are MC data at two different temperatures T = 0.681J//cb and T = 0.682J/£;b; the full curves 
are the best-fit lines. The slopes of the fit lines give v/v = 0.79 ± 0.02 at T = 0.681J/A;b and 
v/u = 0.83 ± 0.02 at T = 0.682J/&B- The error bars from the jackknife method are smaller than 
the symbol sizes. 
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FIG. 11: PRG method to extract the transition temperature T c and the critical exponent v/v of 
the helicity modulus 7. All three curves show j(L, t)/j(L', t) for the same ratio of L/V = 2. The 
x and y coordinates of the intersection of the three curves yield the values of T c and (L / L')~ v / U . 
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FIG. 12: RGT method to extract the transition temperature T c and the critical exponent v of 
the correlation length. The ratio of the slopes of the two fitting lines at T c gives the value of v, 
according to Eq. (|37|) . 
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FIG. 13: Four degenerate eigenmodes for the 2x2x2 unit cell of the ordered state with f = 
(1/2,1/2,1/2), corresponding to the mean-field transition temperature T^ F = v3/2 J/k&. The 
numbers represent the order parameters r\i of the mode, as indicated in Eq. (|39|) . The phases reside 
on the nodes of the lattice. 
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